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At zero chemical potential /i, the order of the temperature-driven quark-hadron transition depends 
on the quark masses m u j and m s . Along a critical line bounding the region of first-order chiral 
transitions in the (m u ^,m s ) plane, this transition is second order. When the chemical potential is 
turned on, this critical line spans a surface, whose curvature at fi = can be determined without 
any sign or overlap problem. Our past measurements on N t = 4 lattices suggest that the region 
of quark masses for which the transition is first order shrinks when ji is turned on, which makes 
a QCD chiral critical point at small jj,/T unlikely. We present results from two complementary 
methods, which can be combined to yield information on higher-order terms. It turns out that the 
ff(H A ) term reinforces the effect of the leading 0(lA 2 ) term, and there is strong evidence that the 
<^(/i 6 ) and ^(/X 8 ) terms do as well. We also report on simulations underway, where the strange 
quark is given its physical mass, and where the lattice spacing is reduced. 
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Figure 1: (Left) Schematic phase transition behavior of Nf = 2 + 1 QCD for different choices of quark 
masses (»z, (i /,ra. s ) at fi = 0. (Middle, Right) Critical surface swept by the chiral critical line as jJ, is turned 
on. Depending on the sign of the curvature c\, a QCD chiral critical point is present or absent [[[J. For heavy 
quarks the curvature has been determined [[|] and the first-order region shrinks with ji . 



1. Introduction 

The fundamental importance of the phase diagram of QCD, as a function of temperature T and 
quark chemical potential ft, makes it the object of several current lattice investigations. It depends 
sensitively on the u, d,s quark masses. At ft = 0, Fig. 1 (Left) summarizes the prevalent under- 
standing of the order of the finite-temperature quark-hadron transition as a function of m u = 
and m s . The physical point lies in the crossover region, separated from the chiral, first-order region 
by a second-order chiral critical line. While the ft = situation is far from settled, it can in princi- 
ple be resolved by manageable increases in computer resources. When fL ^ 0, the complex nature 
of the fermion determinant makes the matter much worse. While finite-/! results, including the 
location of the QCD critical point, have been obtained by reweighting ft = data [^J, assessing the 
reliability of these results is a challenge in itself It appears that the only information that can be 
obtained reliably (i.e. performing thermodynamic and continuum extrapolations) in principle, bar- 
ring an algorithmic breakthrough, is the Taylor expansion of thermodynamic observables in (jl/T) 
about /i = 0. This makes the detection of a finite-/! critical point, characterized by a singularity in 
the free energy, particularly difficult. 

To circumvent this problem, our strategy consists of Taylor-expanding the surface swept by the 
chiral critical line of Fig. 1 (Left). The Taylor expansion of a generic quark mass m c on the chiral 
critical surface, and the associated transition temperature T c , can be written as: 

I^hEl = l+ y au ('^) k (-tL) 2 ', (U ) 

^ = 1 + x> (JLf. (l . 2) 

The sign of ci governs the small-/i behaviour, as illustrated Fig. 1. Our first results [|I|], for the 
Nf = 3 (m s = m Ui d) theory on an 8 x 4 lattice, favored a negative value for c\. In [Q], we presented 
a new numerical method to obtain the c^'s. Here, we combine the two methods and report on our 
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progress towards determining c\ and higher Taylor coefficients (z) on larger lattices; (ii) for the 
Nf = 2+1 theory with physical m s ; (Hi) for the Nf = 3 theory on a finer, N t = 6, lattice. 

2. Extracting the /i-dependence of the critical point 



On the lattice, the Taylor expansion ( |1.2| ) is replaced by that of dimensionless observables: 

f3 c (afn,an) = p c (arn c ,0)+ V cu (am- am c ) k (afL) 21 , (2.1) 

k,i=i 

am c (a\i) = am c Q + £ c' k (afL) 2k . (2.2) 



k=l 



To differentiate between crossover, second- and first-order transitions, we monitor the Binder cu- 
mulant of the quark condensate: 



((8fY) 2 ) : 



B ^ ;LLA 2 > sw=w-(w), (2.3) 



when ((<5yn//) 3 } = 0. On the chiral critical surface, B4 takes value 1.604 as dictated by the 3d Ising 
universality class. It can be expanded as: 

B 4 (am,an) = 1.604+ Y b u (am- am c ) k (ajl) 21 , (2.4) 
k,i=i 

with coefficients satisfying the scaling behaviour bu (L) = fi\L^ k+l ^ v for large L. Having measured 



the first few b^i's by the methods of Sec. 3, we can reconstruct the c^'s eq.(2.2) as: 



j = darrf = dB 4 ( dB 4 \ 1 = ^ 
1 d(a\x) 2 d(a\x) 2 \dam J b\o ' 

c! * = ™ 1u ^2l2 = " T~ ( fc 02 + bi 1 c[ + b 20 c'i) . (2.6) 
2! <i[(c?^) 2 ] 2 />io 



and finally c\ and C2 as: 



, . £1 , 1 dT c (m c (ji),ji) 



N? am c T c (m c ,0) d(n/nT) 
c _ 7l 4 c> 2 K 2 c\ 1 dT c {m c {ii)^) 1 d 2 T c {m c {v),H) 

° 2 N? am c N? am c T c (m c ,0) d(\lj%T) 2 2T c (m c Q ,0) d[(fL / %T) 2 } 2 

3. Two methods to measure B4 derivatives 

B 4 varies steeply with the quark mass, and b\o, b 2 o in eq.(^4j) can be obtained straightforwardly 
from fits of B 4 measured at fL = for different quark masses |JIJ]. Measuring the variation of B4 with 
fL is another matter: B4 is a noisy quantity, its variation is small, and simulating at non-zero (real) fL 
is not feasible. We have used two different, complementary methods to bypass these difficulties [j|] : 
1. We perform simulations at several imaginary values fL = ifLj, where the sign problem is absent, 
and fit our measurements of 54 (jU,-) with a truncated Taylor series in fL 2 . 
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Figure 2: (Left) Comparison of two methods of measuring dB4/d(ajij) 2 on an 8 3 x 4 lattice. The broad 
error band is the fit to imaginary fx data; the data points show the reweighted finite difference quotients, 
obtained with about 4 times fewer statistics. (Right) Finite-size scaling test: data obtained on 8 3 and 12 3 x 4 
lattices show good consistency with the 3d Ising universality class. 



2. We perform simulations at fL = 0, reweight to small values fL = ifLu and measure the finite 
difference quotients kB 4 / A.(afL) 2 , with 

AB 4 dB 4 



lim 



%VoA(fl,li) 2 d(a/i) 2 



(3.1) 

A comparison between the two methods is provided Fig. 2 (Left), on an 8 3 x 4 lattice for Nf = 3. 
The eiTor band is the fit to the finite-/!; data (method 1). The data points are the finite-difference 
quotients (method 2). Consistency between the two methods is observed. The second method is 
clearly more efficient, since the statistics is only 1/4 of the other. This efficiency can be traced to the 
strong cancellation of statistical fluctuations when measuring AB 4 on the ft = and the reweighted 
ensemble. Reweighting itself is done stochastically with a Gaussian-distributed vector Y\ , since the 
reweighting factor is 



p(mi,m 2 ) = "";:; =( e *p . 0.2) 

der> /4 #(t/,jUi) \ v ' 1 V 

Note the small values of {a Hi) 2 in Fig. 2 (Left): they guarantee a good overlap between the fL = 
Monte Carlo ensemble and the reweighted ft = iftj ensemble, and small fluctuations in p. 

Since our 8 3 lattice is not very large (m n L ~ 3.4), we performed a finite-size scaling check by 
comparing with a 12 3 x 4 lattice. Fig. 2 (Right) shows nice consistency with the expected large 
volume universal behaviour, not only for the j-axis intercept yielding Z?oi, but also for the slope 
yielding Z?o2- The result (&02 > like Z?oi) reinforces the finding that the transition weakens and 
turns into a crossover (i.e. B4 increases) as fL is turned on (see eq. (2.4)). 

Finally, we can combine the data from our two methods, since the simulations were performed 
independently and cover different ranges of /»,. A combined fit of the am = 0.0265 data Fig. 3 
shows that (B^afLi) —B^fL = 0))/(afLi) 2 is an alternating series in (afLi) 2 The fit gives 

B 4 {aHi) = B A {fL = 0) - 1.79(14)(aju f ) 2 + 108(27)(^) 4 - 3438(933)(^) 6 + 35954(8876) (a/i,) 8 

(3.3) 
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Figure 3: Combining the two methods: the (a/i,) 2 < 0.01 data come from jj, = reweighting, the (afxi) 2 > 
0.01 from direct ju,- f^ simulations, all at am = 0.0265. Data at larger /I, clearly fall below the G(\lf) 
contribution, indicating a negative /i?-term. The quality of the cubic, S-shape fit favors a positive fxf -term. 
After rotation to real fi, all terms contribute to increasing B4, i.e. pushing the system in the crossover region. 



with a # 2 /d.o.f. of 0.57. The large values of higher-order coefficients indicate that higher-order 
terms become important when /i/r>0.5. However, after rotation to real fL, they all tend to increase 
B4, pushing the system deeper in the crossover region. This only increases the validity of the exotic 
scenario Fig. 1 (Right) up to larger values of fi/T. Conservatively, we trust only the @{fl 2 ) and 
&(fL A ) terms. After continuum conversion following eqs.(2.5-2.8), our final result for Nf = 3 on 
coarse, N t = 4, lattices reads [Q] : 

*M = 1 - 3-3(3) m 2 -47(20) (-tL) 4 - . . . (3.4) 
m c (0) V '\kTJ y 'KkTJ 

4. Towards the Nf = 2+1 continuum limit 



We are currently investigating two reasons why our result eq.(|3.4]) could change qualitatively 
as we consider real QCD. The sign of the curvature could change as we move along the critical line 
away from the degenerate Nf = 3 case. It could also change as we take the continuum limit. 

The first possibility appears unlikely given our current results Fig. 4 (Left), where m s is given 
its physical value on the N t =4 critical line determined in ^ (see Fig. 4 (Right)). Since our pions 
are lighter than in nature, large lattices are required and thereby large computer resources. This is 
achieved, like for the Nj ■ = 3, 8 3 x 4, method 2 case above, by dispatching our simulations over the 
computing Grid. Many independent Monte Carlo runs are performed, all at fL = 0, over a range 
of temperatures near T c , using prioritized scheduling. Current statistics reach 600k thermalized 
configurations. 
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Figure 4: Work in progress: (Left) Nf = 2 + 1 on a 16 x 4 lattice. The simulation point, indicated by the 
leftmost arrow (Right), lies to the left of the physical point, implying that our pions are lighter than in nature. 



The effect of a finer lattice is studied by simulating 18 3 x 6 lattices with Nf = 3 degenerate 
flavors. The current results, Fig. 5 (Left), give opposite signs for &oi using a leading or subleading 
order fit. While the sign of the curvature c\ is consequently not clear, one can already say that \c\ \ 
is not large, £?(20) or less. Thus, the critical surface is almost vertical. 

In addition, another qualitative effect takes place: the ft = critical line, and thereby the whole 
chiral critical surface, moves towards the origin as a — > 0. For instance, the Nf = 3 pion mass on 
the critical line drops from 1.680(4)r c to 0.954(12)7; going from N, = 4 to N, = 6 lattices [§]. 
The first-order region, in physical units, shrinks dramatically as a — ► 0. To compensate this effect 
and maintain a critical point for real QCD at small chemical potentials /x/T < 1, a large positive 
curvature c\ would be needed. We presently do not see it. 

Finally, we note that effective models like PNJL [§] or linear sigma model with simple 
modifications, can reproduce the qualitative features of the chiral critical surface which we observe. 
Nevertheless, let us stress again that our study concerns only the chiral critical surface, swept by 
the jU = chiral critical line as the chemical potential is turned on. Our results do not preclude 
other phase transitions, not connected to the chiral one. 
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Figure 5: Work in progress:(Le/f) Nf = 3 on an 18 3 x 6 lattice. The sign of the intercept boi depends on the 
fitting ansatz. (Middle) As the lattice spacing is reduced, the critical surface moves towards the origin. If a 
critical point exists, its location moves to larger jj,. (Right) If the curvature c\ is negative, higher-order terms 
must increase in magnitude for a critical point to occur. 
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